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On biases in the predictions of stellar population synthesis models 
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ABSTRACT 

Sampling fluctuations in stellar populations give rise to dispersions in observables when a 
small number of sources contribute effectively to the observables. This is the case for a vari- 
ety of linear functions of the spectral energy distribution (SED) in small stellar systems, such 
as galactic and extragalactic Hll regions, dwarf galaxies or stellar clusters. In this paper we 
show that sampling fluctuations also introduce systematic biases and multi-modality in non- 
linear functions of the SED, such as luminosity ratios, magnitudes and colours. In some cases, 
the distribution functions of rational and logarithmic quantities are bimodal, hence complicat- 
ing considerably the interpretation of these quantities in terms of age or evolutionary stages. 
These biases can be only assessed by Monte Carlo simulations. We find that biases are usually 
negligible when the effective number of stars, TV, which contribute to a given observable is 
larger than 10. Bimodal distributions may appear when TV is between 10 and 0. 1 . Predictions 
from any model of stellar population synthesis become extremely unreliable for small TV val- 
ues, providing an operational limit to the applicability of such models for the interpretation of 
integrated properties of stellar systems. In terms of stellar masses, assuming a Salpeter Initial 
Mass Function in the range 0.08 - 120 M Q , TV=10 corresponds to about 10 5 M Q (although 
the exact value depends on the age and the observable). This bias may account, at least in part, 
for claimed variations in the properties of the stellar initial mass function in small systems, 
and arises from the discrete nature of small stellar populations. 

Key words: galaxies: statistics - galaxies: evolution - galaxies: dwarf - galaxies: starburst - 
galaxies: star clusters -methods: statistical 



1 INTRODUCTION AND MOTIVATION 

The comparison of observations to theoretical models is one of the 
basic steps which allows an understanding of natural phenomena. 
The comparison is not always satisfactory and leads to either new 
insights on the nature of the phenomena observed or to revisions 
of the theoretical framework within which the models are made, or 
both. 

In general models have an intrinsic uncertainty, in addition to 
possible systematic effects, and one seeks to minimise both in or- 
der to infer a robust interpretation of the observations. In the case 
of stellar population synthesis, the former uncertainties have often 
been neglected in comparison with the latter, where external errors 
(i.e. uncertai nties on the input assu mpt ions) have been much dis- 
cussed, e.g. ^eitherer et al. (1996) or Bruzual (2001). However, 



the predicted integrated properties of small stellar systems suffer 
from large intrinsic dispersions arising from the very nature of the 
systems: sampling fluctuations from, say, the initial mass function 
(IMF) will give rise to fluctuations in the prop erties and hence to 
an intrinsic dispersion i n the predictions (e.g . fchiosi et al.||l988[ 
Santos & Frogej 1997|; Lan con & Mouhcine ZOOO^ervino et al 



Santos & r/rogel 1199/1; l^an con & Mouhcinel (ZUOOt |Lervino et al. 
ZOOOj , |2001[ |Bruzual||2002| : |cervino et alj |2002| ). This effect has 



been overlooked by most synthesis codes, which usually take an 
analytical IMF and produce predictions which are only valid in the 
limit where the number of stars is formally infinite. 

In addition to this dispersion -which can be quantified and 
properlvtakMMTitojirecam^y^ et 
al. ( |2002| ) - there is also a more subtle effect arising from this in- 
complete sampling. Integrated quantities which scale linearly with 
the number of stars (such as the luminosity in a given photometric 
band, or the overall spectral energy distribution, SED) can be pre- 
dicted for any stellar system from the scaled-down outputs of syn- 
thesis codes with effectively predict properties for very large stellar 
systems, that is, systems which fully sample the distribution func- 
tion of masses. However, for quantities which do not scale linearly 
with the number of stars, such as luminosity ratios, magnitudes, 
colours, equivalent widths, etc, such a scaling cannot be performed 
properly, and the only solution is to make extensive Monte Carlo 
simulations. We can already foresee that biases will arise from the 
incomplete sampling of the underlying distribution function. For 
instance, infrared colours may be dominated, at some ages, by a 
handful of high lu minosity, IR-bright star s. Monte Carlo simula- 
tions calculated by Santos & Frogel ( 1997) show that the average 
colour is a function of the number of stars in the simulation. Simi- 
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Figure 1. Integrated f? — V colour vs. V — K colour of 10 Myr old stellar 
clusters with 1000, 100, 10 and 1 star, respectively, with masses between 2 
and 120 Mq distributed following a Salpeter IMF (triangles). Each triangle 
represents a Monte Carlo simulation with these parameters. The location of 
a stellar population with an infinite number of stars is indicated with the 
symbol. The colours of the sparsely-populated clusters show a multimodal 
distribution in this diagram, and their average values are biased with respect 
to the colours of an infinite stellar population. 



larly, the average optical colours of stellar clusters, at a given age, 
are functions of th e total mass in the cluster, as shown for inst ance 
by|Girardi &Bica| ( |l993| ) and more recently byf 
Bruzuaj ( [200^ 



Girardi 



(200C)and 



This is best illustrated with an example in a simple colour- 
colour diagram. Figure [j] shows the point (indicated by the sym- 
bol) where a 10 Myr-old stellar population with an infinite number 
of stars would be. Each triangle in each panel represents one Monte 
Carlo simulation of a 10 Myr-old stellar population with a total 
number of stars (with masses between 2 and 120 Mq distributed 
with a Salpeter IMF) of 1000, 100, 10 and 1 star, respectively^ 
Figure |l| shows three important effects: 

(a) The distribution of colours in this diagram is clearly bi- 
modal or even multimodal. The points do not cluster around a well- 
defined center, but are very scattered over the entire diagram, at 
odds with the (standard) predicted value for an infinitely-large stel- 
lar population. 

(b) If the average colour is computed for clusters with the same 
number of stars, its value is very different from the one given by 
the standard prediction : the average colour is biased, and the bias 
depends on the total number of stars present in the population. 

(c) The area covered by the points (each representing one 
Monte Carlo simulation) is not a monotonic function of the number 
of stars in the simulation. 

Note that these effects are not limited of course to clusters 
of this age, but do also appear in more evolved populations, since 
these effects are a generic feature of under-sampled clusters^ 

In this paper we provide a first order analytical estimation on 
when and why synthesis models which use an analytical formula- 
tion of the IMF are not be able to reproduce the observed proper- 



1 For reference, jVt o t=1000 stars within 2 and 120 Mq have a weight of 
Aftot ~ 6 X 10 3 Mq, and would the IMF be extended down to 0.08 Mq, 
the mass of the cluster would be M t „+=2.25 xJO 4 Mp 



i . 



2 See for example the analysis by Bruzual ( 2002 ) where the positions of 
LMC clusters are interpreted in terms of Monte Carlo simulations of clus- 
ters of 10 s , 10 4 and 10 3 Mq at different ages. 



ties of stellar systems and therefore when Monte Carlo simulations 
are absolutely required in order to make predictions. We show that 
bias and multimodality features are a generic property of non-linear 
functions of the SED, and will always be present when trying to ap- 
ply population synthesis models to small size stellar populations. 
Their importance cannot be stressed enough: many of the seem- 
ingly 'peculiar' properties of some systems can, at least in part, be 
accounted for by this bias and multimodality effects. 

This is an important issue since odd values of these proper- 
ties often lead to non standard interpretations. For instance, when 
colours or equivalent widths (EWs) of small stellar systems appear 
to be peculiar (that is, peculiar with respect to the values predicted 
by synthesis models of a large number of stars), they are often in- 
terpreted in terms of IMF variability : either the mass cutoffs have 
to be truncated, and/or the slope of the IMF has to be changed in or- 
der to produce the 'correct' colours or EWs. Since the lower mass 
cutoff determines the overall mass-to-light ratio, it is usually the 
upper mass limit which is varied to account for the data. Exam- 
ples of these i nterpretations aboun d in th e recent literature, starting 
perhaps with Shields & Tinsley (1970 who, on the basis of the 
variation of the equivalent width of H/3 in Hll regions across M 
101, inferred a dependence of the upper mass cutoff with metallic- 
ity. Line ratios included in clas sical diagnostic diagrams have also 
been used (e.g. |Plofsson ||l997|, and B resolin, Kennicutt & Garnett 



1999| for the softness parameter) and seem to show that the IMF 



changes in extragjdartir^ftn-ejri^^ & 
Kennicutt |2002| , for an alternative explanation based on sampling 
effects). Colours, and colour-colour diagrams, sometimes includ- 
ing the Ha line, of clusters and starbursts also seem to point to 
variations in the parameters of the IMF, from the LMC (Dolphin 
& Hu nter " |l99% to clusters in the starburst galaxy IC10 (Hunter 
200 lh. 



The common features of all these interpretations which point 
to a variation of the IMF are (i) the systematic use of non-linear 
quantities (EWs, line ratios, colours) and (ii) observations of small 
stellar systems (Hll regions, stellar clusters, starbursts, etc). 

Given the scaling arguments explained above, it is interesting 
to see whether sampling effects could account for this type of ob- 
servations. In this paper we show that whilst the average values 
of the properties (and their dispersions) that scale with the number 
of stars are reasonably well reproduced by the models, non-linear 
functions such as luminosity ratios or colours may be heavily bi- 
ased. It is beyond the scope of this paper to analyse in detail each 
and every observation, but given our results, it appears very likely 
that the 'peculiar' properties observed in some systems can be sim- 
ply understood in terms of the bias and multimodality introduced 
by sampling fluctuations of the IMF in these stellar systems. This 
has profound implications for the interpretation of colours, average 
magnitudes and line ratios in such systems, as well as for the uni- 
versality (or otherwise) of the IMF. These effects seem to be large 
enough that most observations may point to the u niversal i ty of the 
IMF, despite a variety of physical conditions (e.g. Croupa 200^ for 



Mas-Hesse & Kunth 



199S for the case of 



a general review and 
starburst galaxies), and therefore provide a unique insight on star 
formation processes in a variety of environments. 

The structure of the paper is the following: in Sec. §2 we recall 
the scaling properties of synthesis models. In Sec. §3 we analyse 
the statical properties of functions of the SED, showing the biases 
and multimodalities that appear in their distribution functions. We 
discuss these results and their implications in Sec. §4 and §5. 
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2 ON SAMPLING AND SCALING IN POPULATION 
SYNTHESIS 

Population synthesis models often scale predictions to the total 
mass in stars. There are at least two reasons for this. First, a (suit- 
able) fixed mass-to-light ratio can be achieved. Second, the predic- 
tions, which are made for a unit star formation rate (that is, the rate 
at which a given gas mass is transformed into stars), can be scaled 
to other rates, or, equivalently, to other masses in gas or stars. 

The underlying mathematical justification for this scaling is 
simple to see. Take for example the average total luminosity in a 
given band (il =< L >. This is given by the weighted sum of the 
luminosities Li of the JVtot stars present in a population at some 
given age : 



(1) 



where / is the number of initial (zero-age) masses considered such 
that JVtot = X/i=i ni an( ^ ni i s me (integer) number of stars of 
type or class i. For a sample of different size JV' to t, all other pa- 
rameters kept the same (age, metallicity, etc), the following scaling 
applies: 



MAM / JVu 



/j.l(N' tot) / N't 



(2) 



This scaling arises because, at a fixed age, Lj are constant s 
whilst n-i follow a Poisson distribution (see Cervino et al. 2002), 



and the average value of a (weighted) sum of Poisson variables is 
the (weighted) sum of average values. It is this linear scaling that 
allows predictions made for a given JVtot to be scaled up or down 
for a different number of stars. 

However, the scaling is always made with the total mass in 
stars, not the actual number of stars JVtot present in a population. 
Obviously, for a fixed IMF, there is a direct relation between the 
total mass in stars and the number of stars, but there is an important 
difference in that the total mass, Mtot is a random variable when 
JVtot is fixed (and vice versa). 

The PDF of Mtot, as of any function of the IMF, can be ob- 
tained analytically usi ng the characteristic function method (e.g. 
Kendall & Stuart 1977), although the inversion can be cumbersome 
or indeed unworkable. For the purposes of the present paper it is 
enough to show the PDF of Mtot obtained from Monte Carlo sim- 
ulations. 

For JVtot = 1 the PDF of M to t is trivially the IMF itself (see 
Fig. ^bottom right), while for large values of JVtot the central limit 
theorem ensures that it converges to a Gaussian function (Fig. ^| top 
left). For intermediate values of iVtcrt, the distribution has a heavy 
tail at large values of Mtot> as Fig. ^ shows. Note also that this 
tail is not sampled well enough even though 10 4 simulations were 
made, with clusters populated by 1000 stars of masses between 2 
and 120 M© distributed with a Salpeter IMF. The heavy tails are a 
direct result of the poor sampling of the upper mass range, where 
the presence of a few massive stars can be very significant for the 
total mass in a small system. 

The large dispersion in Mtot, at JVtot fixed, shows that synthe- 
sis models which use an analytical 'sampling' of the IMF overlook 
a fraction of the intrinsic dispersion inherent to the nature of the 
problem: small stellar systems will always show a dispersion in 
Mtot event though iVtot is kept fixed, and inversely. So not only 
using JVtot is the natural choice for the scaling, but it also leads 
to take into account one of the factors which produce an intrinsic 
dispersion in the observables. 

The evaluation of the dispersion in Mtot requires to know the 



corresponding PDF, but we can estimate how it scales with iVtot 
from the results of the Monte Carlo simulations. Let us assume that 
the IMF, $(m), is given by a power law: 

dN 
dm 

where A is a normalisation constant. If the IMF is defined between 
some mass cutoffs mi ow and m up , the normalisation constant is 
given by: 



$(m) = Am 



(3) 



.4 



(1 



(4) 



In the case of only one star, the average value of the stellar 
mass is given by 



< m > 



m $(m) dm - 
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(2 - a) (mip c 
and its variance by 



(5) 
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so, square of its relative error is 



}™-L- < m > 2 

1 — a\ 
low / 



(6) 
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< m > 2 
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(3 - a) (mup 1 
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(7) 



For the particular case of a Salpeter IMF, with cutoffs at 2 and 
120 M Q , these values are < m >= 5.897, a(m) 2 = 76.28 and 
a(m)/m = 1.481. 

In general, Mtot is proportional to JVtot- It can be seen that 
the variance, a(M to t) 2 is also proportional to JVtot- Therefore the 
relative error in Mtot is inversely proportional to the square root of 
the number of stars. 

For stellar populations with a large number of stars, the disper- 
sion becomes negligible (e.g., 0.01 per cent for a system of 2x 10 s 
stars) and quantities which depend linearly on the number of stars 
can be normalised safely to Mtot rather than to JVtot (which be- 
comes difficult to measure in such cases!). 



3 DISTRIBUTION FUNCTIONS AND MOMENTS OF 
OBSERVABLES IN SYNTHESIS MODELS 

In order to establish how relevant the effect of sampling is, we have 
followed the evolution of each Monte Carlo simulation of stellar 
clusters of fixed JVtot, and derived their integrated luminosities in 
the V and K bands. 

Using the ensemble of simulations at fixed JVtot, we have ob- 
tained mean values and dispersions for these luminosities, their log, 
the Lv/Lk ratio and the V — K colour. To avoid introducing 
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further parameters in this study, we used the Kurucz (1991 ) atmo- 



sphere mo dels and solar metall icity tracks with standard mass-loss 
rates form Schaller et al. ( 1992 ). For the sake of simplicity, the neb- 



ular contribution to the luminosity has not been taken into account 
in these runs. 

Although Monte Carlo simulations seem the natural, and only , 
way to estimate the PDF of any observable, Cervino et al. (2002) 



showed that it is possible to apply a Poissonian formalism which 
allows the estimation of the first two moments of the PDF analyt- 
ically. This formalism is based on the idea that there is only an 
effective number of stars, J\f, which contribute to a given observ- 
able. Its definition is related to the quasi-Poissonian nature of the 
luminosity PDF. For example, for the observable monochromatic 
luminosity L, its effective number N{L) is given by the expres- 
sion 



1 



ML 



(8) 



Buzzoni 



(1989) 



as first derived by 

Note that this effective number is not an actual number of 
stars, but rather a measure of the weighted number of sources which 
contribute to a given observable. The larger the effective number, 
the better determined will the observable be. Small values of N(L) 
indicate an intrinsically large dispersion, and typically arise when 
a few (real) stars dominate the observable. For instance in young 
clusters a few very massive stars may dominate the overall energy 
budget or the SED. In old stellar populations, a few Asymptotic 
Giant Branch stars (AGB) may change significantly the luminosity. 
Hence in both situations one expects an intrinsically large disper- 
sion in the observable, a direct result of the discrete nature of these 



small populations which sample incompletely the IMF but domi- 
nate the observable. 

Figure ^| shows, for reference, the time evolution of the Ly 
and Lk luminosities (normalised to the number of stars in the clus- 
ter), for a stellar population which fully samples the given IMF, 
that is, with an infinite number of stars in the indicated mass range, 
along with the evolution of their effective number of stars N ', which 
readily allows an estimate of the intrinsic dispersion in these lu- 
minosities via Eq. (^), for any stellar population o f any size and 
any ag e. These were some of the results obtained by 
( |2002| ) 



Cervino et al 



Tt is interesting to see some values of JV. From Fig. H can 
be seen that 5(2.5) x 10 3 stars are needed to obtain J\[Lk) 
(N(Lv)) values larger than 10. This means cluster masses around 
3(1.5) x 10 4 M Q in the mass range 2-120 Mq, and would the IMF 
be extended down to 0.08 Mq, the mass of the cluster would be 
M to t=l. 1(0.6) x 10 5 Mq. 

For young stellar populations, N values for different lu- 
minosities normalised to the stellar cluster mass assuming a 
Salpeter IMF in the mass range 2-120 Mr can b e found in 



http : / / www .laeff.esa.es/users/mcs/SED/ 



In th e case of old ste llar populations N has been computed di- 
rectly by Buzzoni (1989). Additional estimations can be obtained 
frornthe^evaluation^fbrigjita 



zoni (1993 



3 In the case of the work from Worthey (1994) the following relation ap- 
plies: M — M = — 2.51ogA/" — 14.53, where M is the magnitude in a 
given band and M is the brightness fluctuations in magnitudes in the same 
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Figure 4. Comparison between Monte Carlo simulations (Mon) and fully-sampled (Ana) of the average values and dispersions of Ly (left) and Ly (middle). 
The lines correspond to different sizes of stellar populations (iVtot=1000 stars: solid, jVtot=100 stars: dashed, JVtot=10 stars: dot-dashed and iVtot=l star: 
dotted). The top panel on t he rig ht side give the correlation coefficient between both quantities (upper panel) by the simulations (lines) and the estimation 
proposed by Cervino et al (2002) (stars). The ratio between the Monte Carlo and fully-sampled estimates are shown ion the lower panels. 



3.1 Linear functions 

In general, a linear combination of Poisson variables will not be an- 
other Poisson variab le, unless the weights take integer values (see 



Cervino et al 



2002, for a discussion). However, the fact that it is 



a linear combination of random variables allows one to estimate 
exactly its first two moments, even though the remaining moments 
can sometimes be also determined analytically. This is well illus- 
trated by Fig. ^ which shows the comparison of the Monte Carlo 
simulations with the fully-sampled results for the Lv and Lk ob- 
servables. The agreement is very good in all cases, with small devi- 
ations when the number of stars in the cluster decreases. Except in 
rare occasions, the analytical expressions obtained via Eq. are 
within 3 per cent and 5 per cent of the Monte Carlo ones. Not only 
the average values agree at this level (and hence are not biased) but 
also the estimation of the dispersion is very good, over a dynami- 
cal range of three orders of magnitude in the total number of stars. 
Linear functions of the luminosity (such as luminosities integrated 
in a photometric band, or ionising fluxes) are not biased. 

Remarkably, our Poissonian formalism is also able to predict 
the mean value and dispersion for simulations of clusters which 
contain only one star. In this case, the IMF gives us the relative 
weight of the probability Wi to obtain a star with a given luminosity 
U. The mean value, taking each cluster as a one realisation, is 



band, and the term 14.53 = 2.5 log 0.65 X 10 is used to obtain the corre- 
sponding M value normalised to the cluster stellar mass for stars masses in 
the mass range 0.08-1 20 M^ followin g a Salpeter IMF slope, instead the 
normalisation used by |Worthey| ( |l994| ): 10 6 Mq in the mass range 0.1-2 
Mq. 



(9) 



with dispersion 



\L) =\(l 



\2 



2 



(10) 



Eq. (^) implies that when Af is very small ^ifwi >> p? L , 
and so Poisson dispersions and Monte Carlo ones become very sim- 
ilar. 

The highly variable correlation coefficient between Ly and 
Lk is also predicted correctly by the Poisson formalism, as shown 
on Fig. ^j. The agreement is in fact at the 2 per cent level for most 
ages. It shows that a (variable) fraction of stars contribute in lumi- 
nosity to both bands, and therefore the actual dispersions in these 
luminosities must take into acc ount this impo r tant t erm if errors 
are to be evaluated properly. As jCervino et al. ( 2002 ) showed, the 
effect of not including this term is to systematically underestimate 
the dispersion. 

So far so good for the first two moments, but what about the 
actual shape of the PDF of these luminosities? Could the devia- 
tions from the Poisson formalism be accounted for deviations in 
the PDF? Figures ^ and |^ give these PDFs at selected ages in terms 
of the variable x = Lj < L >= L//j,l- Note the logarithmic 
scale on the vertical axis. The vertical dashed lines give the values 

OfX = ^nalytical^MonteCarlo^ as shown above> ^ yery 

close to unity and hence are not biased. It is interesting to note that 
whereas the mean value of the simulations are well reproduced, the 
median value and the most probable (mode) values are shifted to 
the left in all cases. This effect must be taken into account when 
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Figure 5. PDF distributions Q (x) for x = Ly /fi(Ly) at selected ages for a variety of sizes of stellar populations. Note that Af is given only for clusters with 
./Vtot=1000 stars, since it scales linearly with N to t, that is, A/"(A?tot) = Af(10 3 ) X (N to t/I0 3 ) for the other populations. 



the result of standard synthesis models are compared with individ- 
ual clusters. 

At early ages most stars are close to the ZAMS, and hence the 
mass-luminosity relation makes that the PDF in luminosity is very 
similar to the PDF distribution in total mass (compare with Fig. ^). 
As stars evolve, post-MS evolutionary stages are reached, and ob- 
viously the longer the stage, the better sampled it is. This gives rise 
to secondary peaks in the PDF, where stars can accumulate, lead- 
ing to multi-modality in the luminosity distribution. The number 
of such stars will be smaller than 1 per cent of the stellar popula- 
tion of the cluster. When the total number of stars is low, there is a 
non-zero probability that the cluster will have no post-MS. So, for 



a given fixed number of stars, it is possible to have clusters (shar- 
ing the same total number of stars) with or without post-MS stars, 
a clear case of possible bimodality when the observable depends 
on the presence (or otherwise) of post-M S stars. This effect was 
first pointed out in the pioneering work by |Chiosi et"aL ( 1988 ). We 
stress again that the average values are not biased, and that the de- 
viations from the Poisson distribution are relatively small in terms 
of the first two moments only. Obviously higher moments are far 
more sensitive to shorter evolutionary stages, and hence more prone 
to sampling effects. 

We can nevertheless quantify these deviations in terms of J\f. 
A Poisson distribution defined by a parameter smaller than M = 5, 
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Figure 6. Same as Fig. ^ for x = Lk / h(Lk)- 



has a probability larger than 1 per cent to have a number of effective 
sources equal to zero. As mentioned before, the number of effective 
sources is not an actual number of stars, but rather a measure of the 
effect of (under-)sampling in these clusters. For M values between 
5 and 0.1, the probability of having a cluster with zero effective 
sources increases from 1 per cent to 90 per cent. For even smaller 
values of N most clusters will have no effective sources in that 
band, because they do not contain stars in the suitable evolutionary 
phase which contributes to the observable. We therefore expect that 
when M is in the range between about 5 and 0.1 bimodality will be 
most easily observed. 

Again, M is an observable-dependent quantity, so it is possi- 
ble that some clusters have large N values for some quantities but 



small values for others. If these quantities are linear functions of 
the SED or the luminosity, then their average values are not biased. 
Any non-linear function will, on the contrary, be prone to biases 
and may enhance multimodality, as we show in the following sec- 
tions. 

3.2 Logarithmic functions 

In general, standard synthesis codes compute logarithmic quantities 
assuming implicitely that their average value is the logarithm of 
the average. For the logarithm of the luminosity this would be for 
instance 

MlogL =< logL > (11) 
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Figure 7. Comparison of average values and dispersions of log Ly and log Lj- from Monte Carlo simulations and from Eqs. and ( p^ ) (used in synthesis 
models with an infinite number of stars). The line styles and symbols correspond to different sizes of stellar populations (A r tot=10 3 stars : solid lines and 
circles, iVtot =1 2 stars : dashed lines and squares, JVt o t=10 stars : dot-dashed lines and triangles and JVtot=l star : dotted). Symbols correspond to the use of 
Eq. (|18[). Filled symbols correspond to Af(L) values lower than 1 . 




Age [Myr] 

Figure 3. Temporal evolution of the average Ly and Lk luminosities and 
their corresponding effective number of stars, Af, normalised to the actual 
number of stars TVtot ■ This is the evolution predicted for an infinite popula- 
tion of stars, fully sampling a Salpeter IMF between 2-120 M©. 



and the corresponding dispersion would be evaluated as 

(loge) 2 



o- 2 (logL) = (loge) 2 -j 



M{L) 



(12) 



This is obviously not correct in general, since the assumption 
relies implicitely on an infinitely narrow-peaked PDF for L, which 
is seldom the case. 

Let us assume a sample of observations of a definite-positive 
random variable L described by a PDF V(L). The expectation 
value [ic of any function of this variable C = C(L) is 



fic = < C > = 



C(L)V(L)dL. 



(13) 



Take for instance the case where L is the luminosity in some 
photometric band. The corresponding magnitude is M = a + 
bin L, with a some zero point and b = — 2.5 log 10 e. Its average 
value is 



fj, M = < M > = 



'0 

a + b 



MP(L)dL 

(In L)V(L) dL. 



(14) 



With the new variable x = L/fiL, whose PDF is Q(x) such 
that of course J Q(x)dx — J V(L)dL — 1, this results in 



fJ-M = a + b In )1l + b 



(\nx) G{x) dx. 



(15) 



Therefore, unless G{x) — S(x — 1) (where 8(x) is Dirac's 
Delta function) -that is, when the PDF of L is defined only at hl 
- the average of the logarithmic quantity will not be the logarithm 
of the average value of the variable. The proper computation of 
the Q(x) functions can be only performed with a massive number 
of Monte Carlo simulations. From this properly-sampled Q(x) one 
could extract samples of x for any sample of size iVtot, but this 
is obviously unpractical. If this is nevertheless not done, then the 
estimated values are very likely to be biased. 

A rough estimate of the bias can be evaluated by expanding in 
Taylor series. For the magnitude M this is 



M = a + b In fi L + b 



L 



2 Hi 



(16) 
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Figure 8. PDF of log L y at selected ages. The vertical dashed line gives the position of log < L y > and hence of the bias. Note the variable amount of bias 
as a function both of age and size of the stellar population. 



so that, to first order, the mean value of the magnitude, and its vari- 
ance are 



N{L) \ 2JV(L) 



+ ■ 



(18) 



2 

fiM = a + b In ^ — b —kr- + 

2/4 



= a + b In hl — 



2Af{L) 



+ 



(17) 



2 



4 m! 



To first order, the bias is then 2 a/(l) > so ^ N(L) is large 
enough, Eqs. ( |ll| ) and ( |l2[ ) are good approximations. In many cases 
N(L) will be small and the bias important. Eq. jl8| ) gives a lower 
limit to the bias (since higher order corrections were not taken into 
account) and the only way to estimate the bias is, inevitably, via 
Monte Carlo simulations. This is done in Fig. ^ where this analyt- 
ical estimate is compared with the Monte Carlo ones (middle and 
right panels). 

Alternatively, in terms of difference in magnitude 
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AM L = n(M L ) -Ml = 2.5(logL - WogL ), (19) 

we see that the bias can reach some 0.5 magnitudes in My and 
2.5 magnitudes in Mk for clusters with 100 stars between 2 and 
120 Mq, and 2 magnitudes in My and ~ 5 magnitudes in Mk 
for groups of 10 stars in that mass range. We can also see that us- 
ing Eqs. ( p"l| ) and ^ give systematically wrong results, because 
they predict systematically too large values. This is easy to under- 
stand since they assumed a fully-sampled IMF, were massive stars 
are (comparatively) over represented with respect to Monte Carlo 
samples of small size. The estimates given by Eq. ( |l8| ) are indi- 
cated with symbols in Fig. R Although they fail to give the correct 



correction (due to the lack of higher order terms) in the bias, they 
provide a rough indication of the underestimation in the dispersion. 

Figs. |I] and ^| present the PDFs of log Ly and log Lk respec- 
tively. In clusters with AT(Ly) or N(Lk) larger than 10 (see Figs. 
@and|), the PDFs are unimodal, as expected. However, for smaller 
values of N the situation is entirely different. In this case, the clus- 
ter is dominated by main sequence stars, which do not sample prop- 
erly the IMF, given their small numbers, and in particular the upper 
mass end. 

At 5.5 and 10 Myr, the Ly distribution of clusters with 100 
stars have Af(Ly) values of 0.7 and 2.5 respectively (see Fig. 
The associated PDFs of log Ly show a bimodal structure with two 
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peaks of similar amplitude (Fig. The effect is more extreme in 
the case of log Lk (Fig. ^) were bimodality is present even for 
clusters with 1000 stars (in the 2-120 Mq mass range) at 5.5 Myr 
with JV(L K ) = 2.3 (see also the PDFs at 5.5 and 10 Myr for 
clusters with 100 stars, with N(Lk) values of 0.23 and 0.54 re- 
spectively). Bimodality with peaks of roughly similar amplitude 
implies that real clusters should have (log)luminosities peaking 
around (roughly) two values (barring selection effects, etc). One 
could naively think that these are two distinct populations of clus- 
ters where in fact there is only one. Since the total number of stars 
is the same, one could be driven to conclude, mistakenly, that the 
IMF in these two 'populations' is different, while in fact it is just a 
sampling effect. 

A final comment on simulations of clusters with one star, be- 
cause their PDFs help to understand the origin of multimodality and 
bias. At 0.5 Myr the PDF is roughly a power law (the convolution 
of a IMF power law with the power-law relation between mass and 
luminosity for main sequence stars). Figs. |s| and ^ indeed show a 
linear relation in these log-log plots. As time goes by, the more lu- 
minous stars (i.e. the most massive ones) leave the main sequence, 
and become more luminous than the brighter main sequence stars. 
The PDFs show two distinct regimes: (i) a power-law for the less 
massive (and thus dimmer) stars, and (ii) secondary peaks for the 
longer post-MS evolutionary phases. When more massive clusters 
are considered, these post-MS peaks grow in importance, become 
better sampled, and eventually the bias disappears. 

3.3 Rational and log-rational functions 

Rational functions of luminosities are even more important in the 
interpretation of the properties of small stellar systems because, 
naively, they seem to be independent of normalisation problems, 
since the usual scaling with mass or number of stars cancels out. 
But this may not be necessarily the case, and, in addition, their dis- 
persion due to the uneven sampling of the IMF can be dramatically 
large. 

In standard synthesis models it is assumed that for a rational 
function u — x/y its average value is given by 



flu — fix x Pi/y — fin I fly-, 
and its dispersion is given by: 



°t , a v o cov(a,y) 

2 ' 2 
fix fly fix fly 



(20) 



1 1 

+ 



- 2 



Af(a:) N{y) ^N{x)N{y) ^( U Y 



(21) 



cov (x,y) 



is the correlation coefficient between x 



where p(x,y) 
and y. 

Again, this would be correct if the distribution of u would be 
a very narrow peak, but in the general case this is not correct. For a 
function R — 1/L, the average value is given by 



flR 



-V(L)dL= — 
L fih 



io ^ Jo 
where x = L / fiL- Again, unless G{x) 

fiR = fll/L # 1/fiL, 
and, more generally, 
H(L\/La) fi Ll /fi L2 . 



■ G(x) dx, 
= 5(x - 1) 



(22) 



(23) 



(24) 



With a second-order Taylor expansion one gets, to first order 



flR = 



2 

<?R 



1 4 

flL \ fit 

2 / „2 



— ii 

flL 



jV(L) 



fi\M{L) \ M{L) 



and, in general, for a rational function u — x/y 



(25) 
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p(x,y) 



M(v) y/M[x)M[y) 
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fix (Ty cov(x, y) x 
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i4 
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+ ... 

p(x,y) 



■A/"(as) A%) yJN{x)N{y) 
( _J_ p{x,y) \ + 



(26) 



To first order approximation, the bias is then of the order of 



i 



P(x,y) 



— ^== i , and the correct estimation of the mean 

y/M(x)N{y)) 

value of a ratio requires the covariance between the numerator and 
denominator. 

Again, we illustrate here these effects with the simple ratio 
Lv I Lk in Fig. [l^, where Eqs. ( ^o| ) and ( pl[ ) produce a very se- 
vere bias, even larger than one order of magnitude, both in the 
average and in the dispersion. The symbols show the results us- 
ing Eq. ( [jo] ) and they have the same meaning as in Fig. ^ except 
for filled symbols, which now correspond to M{Lv) and M(Lk) 
values lower than 10. In this case, the use of the Eq. ( po] ) only pro- 
duces a marginally better result than Eq. (po[). Note that in the case 
of a(Ly I Lk) there are fewer points computed with Eq. (^) than 
in the case of the mean value. This is just because, for some val- 
ues, Eq. (^) produces negative a 2 (Lv / Lk) values and the use of 
these equations is not valid. It is interesting to see that the bias can 
be positive or negative, a direct consequence of the time variability 
of the correlation coefficient between the two luminosities. 

The corresponding PDFs are given in Fig. 111]. In this case, 
bimodality is apparent for a wide range of situations, and again, it 
is due to a secondary peak with larger values of the Lv / Lk ratio. 
The bias may increase or decrease the average value, as shown 
by the vertical dashed lines. It is remarkable to see that the actual 
dispersion for the more massive clusters, at 10 Myr for example, 
is actually larger than the dispersion in less massive clusters. 
Again this is due to the presence of very massive stars, better 
sampled in these massive clusters, which have in comparison a 
much smaller flux in the K band, and hence a very large Lv / Lk 
ratio. The behaviour of the bias with the number of stars is clearly 
non-monotonic, as shown also in Fig. [l^. 
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Figure 10. Mean values and dispersions of Ly LpK to Monte Carlo simulations and in fully-sampled synthesis (via Eqs. ( p(i[ ) and (pT|), as used in standard 
synthesis models). Symbols and lines as in Fig. M except for filled symbols, that now correspond to J\f(Ly) and J\f(Ljc) values lower than 10. Note that the 
a value for clusters with 10 3 stars is larger at 5 Myr than the a value for 'clusters' with one star. It is due to the PDF (Fig. |ll|) of the clusters with one star, see 
text for explanations. 



Finally, in the case of log-rational functions (such as integrated 
colours) there will be a priori a double bias, arising from the use of 
two non-linear functions. And yet most codes will still assume that 
the colours of a stellar population will not depend on its total mass. 
We take again as an example the V — K colour. To first order, 



the predictions of synthesis models must be corrected toward blue 
colours in order to reproduce the mean value of a set of observed 
clusters. At older ages, the bias may turn to bluer or redder values 
depending on the age, that is, on the proportion of post-MS versus 
main sequence stars. 



Hv-k = 



2 

°V-K 



b In 



1 



M{L V ) N{L K ) 



+ . 



+ 



- 2- 



Af(Lv) N{L K ) y/M(L v )N{L K 



N{L V ) M{L K ) 



+ . 



(27) 



Note that in this case the bias does not dependent on 
the correlation coefficient, and it has a value of approximately 

I ( TnjTy) ~ a/(lk) ) ' S mce redder wavelengths have a lower M 
valu e (in this range of ages) than shorter wavelengths (see Cervino 
et al. [2002| for details), this approximation for the bias indicates 
that predictions of synthesis models must be corrected with a blue 
shift in order to get back to the actual mean value of a set of ob- 
served clusters. Note that for older clusters (say, above 1 Gyr or so), 
the trend is reversed since redder wavelengths (cooler evolutionary 
phases) become better sampled. 

The situation is illustrated in Figs. [l2| and [l3|. The top right 
panel of Fig. [l2| shows the mean V — K colour for the different sets 
of simulations. In this case, the difference between actual colours 
and those predicted by fully-sampled synthesis models can reach 
more than 3 magnitudes, even for the most massive clusters. 

In the case of young star forming regions (younger than some 
20 Myr), red colour s have lower N values than blue ones (see 
Cervino et al. 2002, for details). The resulting bias shows that 



Another interesting result shown in Figs. 0and0 is the cor- 
responding a value. In the case of Gaussian distributions, a cor- 
responds to the half width of maximum at about 60 per cent of 
the full height. In our case, the PDFs are not Gaussians, but a still 
gives us the width of the distribution at some (undefined) height. 
The Monte Carlo simulations show that at some ages (when post- 
MS stars appear in the clusters), the a value is maximum for the 
simulations with 100 stars and decreases for even fewer stars. This 
value is also related with the area covered by the simulations in 
Fig. [j], where the 'dispersion' in the diagram is maximum for sim- 
ulations with 100 stars. As we have already seen, this maximum 
a value is also related with a larger probability of obtaining bi- 
modal distributions. The reason of this behaviour is that most of 
the clusters with a smaller number of stars will not have the post- 
MS phases sampled at all. This is also consistent with the lack of 
evolution in the Ly / Lk ratio and the V — K colour. So, the proba- 
bility of finding a cluster with post-MS falls below the width of the 
level defined by a. For more populated clusters, almost all of them 
will be reproduced correctly by standard synthesis models with a 
non-zero number of post-MS stars. 



It is worth emphasizi ng that th i s beh aviour is also present in 
simulations performed by Bruzual (2002). In his diagrams there 



is a maximum dispersion for the simulations with 10 M@ that 
decreases for simulations with 10 3 Mq. Also, the bias effect is 
clearly present in his Monte Carlo simulations since the colours of 
these simulations do not agree with his analytical predictions. 
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4 DISCUSSION AND IMPLICATIONS 

The comparison of observations with the results of a population 
synthesis method is in fact an inverse problem : Given the data, 
what are the models which best fit the data, given the errors in the 
data? It is seldom the case that errors within the models are also 
taken into account. It is also important to analyse whether the best 
fitting models are also good fits, and whether the solution is unique 
or not. In the case of synthesis models this procedure is often ne- 
glected and a variation in the input parameters is much preferred. 

We have shown in the previous sections two major problems. 
First, the usual scaling with the total mass of the stellar popula- 
tion should not be made when dealing with small systems (dwarf 



galaxies, starbursts, stellar clusters, pixels, etc). Second, observ- 
ables which are non-linear functions of the SED (or equivalently of 
the monochromatic luminosity) may present not only multimodal 
distribution functions, but also a systematic bias in the average val- 
ues of these observables, with respect to the values predicted by 
codes which fully sample the IMF. 

These effects have been overlooked in the past, because most 
synthesis codes do not provide an estimate of the intrinsic disper- 
sion of the observables. Whilst this was easy to understand when 
dealing with the integrated properties of massive stellar popula- 
tions, such as spiral and elliptical galaxies, this becomes a serious 
problem if such codes are applied to smaller stellar systems. For- 
mally, the first moment of the PDF of a given observable cannot 
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Figure 12. Mean values and dispersions of the integrated V — K colour from Monte Carlo simulations, and their comparison to simulation results obtained 
directly from the corresponding band luminosities. Symbols as in Fig. llOl 



be scaled down properly. A possible solution would be to have the 
PDF of the observable and then extract samples of small size from 
that PDF. In this way the statistical properties of the observable, 
in small samples, can be analysed and compared with the observa- 
tions. 

The only proper solution to deal with the two problems we 
have discussed in this paper, multimodality and bias, is therefore 
to perform Monte Carlo simulations to construct the PDF of the 
observable. This is often unpractical, and we have shown that the 
quasi-Poisson formalism that we have developed can account for 
the statistical properties of the first two moments at a very good 
approximation level. Although the formalism cannot account for 
higher order moments, and thus to allow for strong deviations in 
the PDFs, it does not suffer from systematics and can easily (and 
should) be included in any synthesis code. 

As we mentioned before, it has often been the case that 
the easiest solution to account for 'peculiar' observations was to 
change some of the input parameters of the synthesis code, and in 
particular the properties of the IMF. It is well beyond the scope of 
this paper to analyse in detail each and every claim for variations 
of the IMF in small stellar systems, but we will illustrate here how 
sampling effects may have been mistakenly overlooked. We take 
again the V — K colour, a non-linear function of the SED, as an ex- 
ample. Figure |l4| (right panel) shows its evolution with time when 
the upper mass limit of the IMF is changed. It is readily apparent 
that if an observation gives a very blue V — K colour, and all other 
parameters are the same (such as mass, luminosity, etc) within a 
small range, then it is tempting to hastily conclude that the upper 
mass limit in that particular system has to decrease. But in fact, as 
Fig. [l2| shows, not only the average colour changes with the mass 
of the system, but also a large dispersion is expected, arising only 



from sampling fluctuations. Similarly, Fig. (left panel) shows the 
effect of varying the slope of the IMF. Again, a too blue colour may 
be interpreted as evidence for a steeper slope. While our results can- 
not confirm or reject possible IMF variations in these systems, they 
however clearly indicate that part (if not all) of the dispersion in the 
V — K colour may be accounted for by the incomplete sampling 
of a single, universal IMF. 



Two final comments. First, in some cases there is no need for a 
bias to be present to explain some peculiar properties of small stel- 
lar systems. All the values quoted in Fig. Il4 lie in fact in the pre- 
dicted band quoted by the maximum and minimum V — K values 
obtained from the simulations of clusters with 1000 stars (where 
the bias is small). The very large dispersion, an expected intrin- 
sic property due to incomplete sampling, may account for many of 
the observed values, with no need to invoke special properties. As 
an example, in high-metallicity regions with Wolf-Rayet stars, the 
dispersion of synthesis models around thei r mean value is enough 



t o exp lain the observations, as shown by Bresolin & Kennicutt 
(§002). 



Second, another interesting effect concerns the ages of these 
systems, as inferred from non-linear functions of the luminosity. 
Because these observables may have a wide, multimodal distribu- 
tion function, the intrinsic dispersion which may be observed in 
a given sample is not necessarily a reflection of a variation in the 
properties of individual members of that sample. It appears very 
dangerous to infer ages or mass variations from these small sys- 
tems. 
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5 CONCLUSIONS 

In this work we have found important effects which have been over- 
looked in the application of population synthesis codes to small 
stellar systems such as clusters, starbursts, dwarf galaxies or pix- 
els. Most current synthesis codes predict the average values of ob- 
servables for an infinitely large population of stars, which samples 
perfectly a given IMF. Extrapolating these predictions to small stel- 
lar systems may be misleading. 

First, scaling observables with the mass in stars is not correct 
for small systems. The proper scaling should be performed with the 
total number of stars present in the system. 

Second, the distribution functions of observables which are 



non-linear functions of the monochromatic luminosity (or SED) 
may show multimodality. This implies that for a given population, 
of fixed size, and all other parameters kept the same, there will be 
a large dispersion in the observable, whose average value will not 
necessarily be the same as the one predicted when the IMF is fully 
sampled. Most synthesis codes will therefore lead to biased predic- 
tions when applied to small stellar systems. 

Third, we have also shown that bimodality is a natural effect 
that appears when stellar systems are affected by sampling (i.e. the 
number of stars in the clusters is not enough to sample completely 
the IMF). It is related to the evolutionary speed of different masses, 
or in more graphical terms, to the absence/presence of stars in some 
evolutionary phases. This effect was already shown 14 years ago by 
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Figure 14. Evolution of the V — K colour for different IMF slopes (left, assuming that M up = 120 Mq) and M up values (right, assuming a Salpeter IMF 
slope). The areas in dark gray shade corresponds to the 90 per cent confidence interval for clusters with 1000 stars. The areas in light gray shade corresponds 
to the maximum and minimum values obtained in simulations of clusters with 1000 stars. 



Chiosi et al. ( 1988 ). We have now established the range of N values 
(from ss 5 to 1) where bimodality appears. An additional result is 
that bimodality will be more easily observed when infrared colours 
are used since such colours are more affected by the presence of 
short-lived post-MS phases. 

Fourth, while we cannot confirm or reject the hypothesis of 
possible variations in the properties of the IMF in small stellar sys- 
tems, these results suggest that an important fraction of the disper- 
sion observed in non-linear observables is not necessarily evidence 
for IMF variations, but can rather be accounted for by sampling 
fluctuations. 

Finally, we have also derived an operational limit to the valid- 
ity of standard synthesis codes. There is a critical effective number 
of stars (unrelated to an actual number of stars), with a value of 
around 10, which defines a boundary where multimodality may ap- 
pear, i.e. cluster masses larger than 10 5 Mq for a Salpeter IMF 
in the range 0.08 - 120 Mq (but the exact value for the cluster 
masses depends on the age and the observable). Below this critical 
value, the results of the codes must be taken with caution, since not 
only they can be biased, but also they may underestimate the actual 
dispersion of the observable, and Monte Carlo simulations are re- 
quired. The reason for the existence of this critical value is directly 
related to the presence (or otherwise) of stellar evolutionary stages 
which are not fully sampled in small systems. 

We encourage the inclusion of the calculation of this effective 
number of sources for a given observable in all synthesis codes, so 
that observers (and theoreticians alike) can estimate the expected 
dispersion (and possible biases) in observables. 
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